Dada la siguiente ecuación diferencial:
$$ \dot{x} = -x $$queremos obtener la respuesta del sistema que representa, es decir, los valores que va tomando $x$.
Si analizamos esta ecuación diferencial, podremos notar que la solución de este sistema es una función $\varphi(t)$, tal que cuando la derivemos obtengamos el negativo de esta misma función, es decir:
$$ \frac{d}{dt} \varphi(t) = -\varphi(t) $$y despues de un poco de pensar, podemos darnos cuenta de que la función que queremos es:
$$ \varphi(t) = e^{-t} $$Sin embargo muchas veces no tendremos funciones tan sencillas (ciertamente no es el caso en la robótica, donde usualmente tenemos ecuaciones diferenciales no lineales de orden $n$), por lo que en esta práctica veremos algunas estrategias para obtener soluciones a esta ecuación diferencial, tanto numéricas como simbolicas.
El método de Euler para obtener el comportamiento de una ecuación diferencial, se basa en la intuición básica de la derivada; digamos que tenemos una ecuación diferencial general:
$$ \frac{dy}{dx} = y' = F(x, y) $$en donde $F(x, y)$ puede ser cualquier función que depende de $x$ y/o de $y$, entonces podemos dividir en pedazos el comportamiento de la gráfica de tal manera que solo calculemos un pequeño pedazo cada vez, aproximando el comportamiento de la ecuación diferencial, con el de una recta, cuya pendiente será la derivada:
«Método de Euler» por Vero.delgado - Trabajo propio. Disponible bajo la licencia CC BY-SA 3.0 vía Wikimedia Commons.
Esta recta que aproxima a la ecuación diferencial, podemos recordar que tiene una estructura:
$$ y = b + mx $$por lo que si sustituimos en $m$ la derivada y $b$ con el valor anterior de la ecuación diferencial, obtendremos algo como:
$$ \overbrace{y_{i+1}}^{\text{nuevo valor de }y} = \overbrace{y_i}^{\text{viejo valor de }y} + \overbrace{\frac{dy}{dx}}^{\text{pendiente}} \overbrace{\Delta x}^{\text{distancia en }x} $$pero conocemos el valor de $\frac{dy}{dx}$, es nuestra ecuación diferencial; por lo que podemos escribir esto como:
$$ y_{i+1} = y_i + F(x_i, y_i) \Delta x $$Resolvamos algunas iteraciones de nuestro sistema; empecemos haciendo 10 iteraciones a lo largo de 10 segundos, con condiciones iniciales $x(0) = 1$, eso quiere decir que:
$$ \begin{align} \Delta t &= 1 \\ x(0) &= 1 \\ \dot{x}(0) &= 1 \end{align} $$
In [ ]:
x0 = 1
Δt = 1
# Para escribir simbolos griegos como Δ, tan solo tienes que escribir su nombre
# precedido de una diagonal (\Delta) y teclear tabulador una vez
F = lambda x : -x
x1 = x0 + F(x0)*Δt
x1
In [ ]:
x2 = x1 + F(x1)*Δt
x2
Crea codigo para una iteración mas con estos mismos parametros y despliega el resultado.
In [ ]:
x3 = # Escribe el codigo de tus calculos aqui
In [ ]:
from pruebas_2 import prueba_2_1
prueba_2_1(x0, x1, x2, x3, _)
Momento... que esta pasando? Resulta que este $\Delta t$ es demasiado grande, intentemos con 20 iteraciones:
$$ \begin{align} \Delta t &= 0.5 \\ x(0) &= 1 \end{align} $$
In [ ]:
x0 = 1
n = 20
Δt = 10/n
F = lambda x : -x
x1 = x0 + F(x0)*Δt
x1
In [ ]:
x2 = x1 + F(x1)*Δt
x2
In [ ]:
x3 = x2 + F(x2)*Δt
x3
Esto va a ser tardado, mejor digamosle a Python que es lo que tenemos que hacer, y que no nos moleste hasta que acabe, podemos usar un ciclo for y una lista para guardar todos los valores de la trayectoria:
In [ ]:
xs = [x0]
for t in range(20):
xs.append(xs[-1] + F(xs[-1])*Δt)
xs
Ahora que tenemos estos valores, podemos graficar el comportamiento de este sistema, primero importamos la libreria matplotlib:
In [ ]:
%matplotlib inline
from matplotlib.pyplot import plot
Mandamos a llamar la función plot:
In [ ]:
plot(xs);
Sin embargo debido a que el periodo de integración que utilizamos es demasiado grande, la solución es bastante inexacta, podemos verlo al graficar contra la que sabemos es la solución de nuestro problema:
In [ ]:
from numpy import linspace, exp
In [ ]:
ts = linspace(0, 10, 20)
plot(xs)
plot(exp(-ts));
Si ahora utilizamos un numero de pedazos muy grande, podemos mejorar nuestra aproximación:
In [ ]:
xs = [x0]
n = 100
Δt = 10/n
for t in range(100):
xs.append(xs[-1] + F(xs[-1])*Δt)
ts = linspace(0, 10, 100)
In [ ]:
plot(xs)
plot(exp(-ts));
Este método funciona tan bien, que ya viene programado dentro de la libreria scipy, por lo que solo tenemos que importar esta librería para utilizar este método.
Sin embargo debemos de tener cuidado al declarar la función $F(x, t)$. El primer argumento de la función se debe de referir al estado de la función, es decir $x$, y el segundo debe de ser la variable independiente, en nuestro caso el tiempo.
In [ ]:
from scipy.integrate import odeint
In [ ]:
F = lambda x, t : -x
In [ ]:
x0 = 1
ts = linspace(0, 10, 100)
xs = odeint(func=F, y0=x0, t=ts)
In [ ]:
plot(ts, xs);
Grafica el comportamiento de la siguiente ecuación diferencial.
$$ \dot{x} = x^2 - 5 x + \frac{1}{2} \sin{x} - 2 $$Nota: Asegurate de impotar todas las librerias que puedas necesitar
In [ ]:
ts = # Escribe aqui el codigo que genera un arreglo de puntos equidistantes (linspace)
x0 = # Escribe el valor de la condicion inicial
# Importa las funciones de librerias que necesites aqui
G = lambda x, t: # Escribe aqui el codigo que describe los calculos que debe hacer la funcion
xs = # Escribe aqui el comando necesario para simular la ecuación diferencial
plot(ts, xs);
In [ ]:
from pruebas_2 import prueba_2_2
prueba_2_2(ts, xs)
Y por ultimo, hay veces en las que incluso podemos obtener una solución analítica de una ecuación diferencial, siempre y cuando cumpla ciertas condiciones de simplicidad.
In [ ]:
from sympy import var, Function, dsolve
from sympy.physics.mechanics import mlatex, mechanics_printing
mechanics_printing()
In [ ]:
var("t")
In [ ]:
x = Function("x")(t)
x, x.diff(t)
In [ ]:
solucion = dsolve(x.diff(t) + x, x)
solucion
Implementa el codigo necesario para obtener la solución analítica de la siguiente ecuación diferencial:
$$ \dot{x} = x^2 - 5x $$
In [ ]:
# Declara la variable independiente de la ecuación diferencial
var("")
# Declara la variable dependiente de la ecuación diferencial
= Function("")()
# Escribe la ecuación diferencial con el formato necesario (Ecuacion = 0)
# adentro de la función dsolve
sol = dsolve()
sol
In [ ]:
from pruebas_2 import prueba_2_3
prueba_2_3(sol)
Si ahora queremos obtener el comportamiento de una ecuacion diferencial de orden superior, como:
$$ \ddot{x} = -\dot{x} - x + 1 $$Tenemos que convertirla en una ecuación diferencial de primer orden para poder resolverla numericamente, por lo que necesitaremos convertirla en una ecuación diferencial matricial, por lo que empezamos escribiendola junto con la identidad $\dot{x} = \dot{x}$ en un sistema de ecuaciones:
$$ \begin{align} \dot{x} &= \dot{x} \\ \ddot{x} &= -\dot{x} - x + 1 \end{align} $$Si extraemos el operador derivada del lado izquierda, tenemos:
$$ $$\begin{align} \frac{d}{dt} x &= \dot{x} \\ \frac{d}{dt} \dot{x} &= -\dot{x} - x + 1 \end{align}$$ $$O bien, de manera matricial:
$$ \frac{d}{dt} \begin{pmatrix} x \\ \dot{x} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & -1 \end{pmatrix} \begin{pmatrix} x \\ \dot{x} \end{pmatrix} + \begin{pmatrix} 0 \\ 1 \end{pmatrix} $$Esta ecuación ya no es de segundo orden, es de hecho, de primer orden, sin embargo nuestra variable ha crecido a ser un vector de estados, por el momento le llamaremos $X$, asi pues, lo podemos escribir como:
$$ \frac{d}{dt} X = A X + B $$en donde:
$$ A = \begin{pmatrix} 0 & 1 \\ -1 & -1 \end{pmatrix} \quad \text{y} \quad B = \begin{pmatrix} 0 \\ 1 \end{pmatrix} $$y de manera similar, declarar una función para dar a odeint.
In [ ]:
from numpy import matrix, array
In [ ]:
def F(X, t):
A = matrix([[0, 1], [-1, -1]])
B = matrix([[0], [1]])
return array((A*matrix(X).T + B).T).tolist()[0]
In [ ]:
ts = linspace(0, 10, 100)
xs = odeint(func=F, y0=[0, 0], t=ts)
In [ ]:
plot(xs);
Implementa la solución de la siguiente ecuación diferencial, por medio de un modelo en representación de espacio de estados:
$$ \ddot{x} = -8\dot{x} - 15x + 1 $$Nota: Tomalo con calma y paso a paso
- Empieza anotando la ecuación diferencial en tu cuaderno, junto a la misma identidad del ejemplo
- Extrae la derivada del lado izquierdo, para que obtengas el estado de tu sistema
- Extrae las matrices A y B que corresponden a este sistema
- Escribe el codigo necesario para representar estas matrices
In [ ]:
def G(X, t):
A = # Escribe aqui el codigo para la matriz A
B = # Escribe aqui el codigo para el vector B
return array((A*matrix(X).T + B).T).tolist()[0]
ts = linspace(0, 10, 100)
xs = odeint(func=G, y0=[0, 0], t=ts)
In [ ]:
plot(xs);
In [ ]:
from pruebas_2 import prueba_2_4
prueba_2_4(xs)
Sin embargo, no es la manera mas facil de obtener la solución, tambien podemos aplicar una transformada de Laplace, y aplicar las funciones de la libreria de control para simular la función de transferencia de esta ecuación; al aplicar la transformada de Laplace, obtendremos:
$$ G(s) = \frac{1}{s^2 + s + 1} $$
In [ ]:
from control import tf, step
In [ ]:
F = tf([0, 0, 1], [1, 1, 1])
In [ ]:
xs, ts = step(F)
In [ ]:
plot(ts, xs);
Modela matematicamente la ecuación diferencial del ejercicio anterior, usando una representación de función de transferencia.
Nota: De nuevo, no desesperes, escribe tu ecuación diferencial y aplica la transformada de Laplaca tal como te enseñaron tus abuelos hace tantos años...
In [ ]:
G = tf([], []) # Escribe los coeficientes de la función de transferencia
xs, ts = step(G)
plot(ts, xs);
In [ ]:
from pruebas_2 import prueba_2_5
prueba_2_5(ts, xs)
Modela matematicamente la suspensión de un automovil de masa $m = 1200 kg$, considera que los resortes de su suspensión tienen una constante $k = 15,000 \frac{N}{m}$ y un amortiguador con constante $c = 1,500 \frac{N s}{m}$.
Gráfica el comportamiento del sistema ante una fuerza $F = 1 N$.